home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / DDJMAG / DDJ8609.ZIP / ASH.SEP next >
Text File  |  1986-09-30  |  18KB  |  664 lines
  1. /* SPLINE.C - Interpolate Smooth Curve
  2.  *
  3.  * Version 2.00        December 25th, 1985
  4.  *
  5.  * Modifications:
  6.  *
  7.  *   V1.00 (85/11/01)   - beta test release
  8.  *   V2.00 (85/12/25)   - general revision
  9.  *
  10.  * Copyright 1985:      Ian Ashdown
  11.  *                      byHeart Software
  12.  *                      620 Ballantree Road
  13.  *                      West Vancouver, B.C.
  14.  *                      Canada V7S 1W3
  15.  *
  16.  * This program may be copied for personal, non-commercial use
  17.  * only, provided that the above copyright notice is included in
  18.  * all copies of the source code. Copying for any other use
  19.  * without previously obtaining the written permission of the
  20.  * author is prohibited.
  21.  *
  22.  * Synopsis:    SPLINE [option] ...
  23.  *
  24.  * Description: SPLINE takes pairs of numbers from the standard
  25.  *              input as abscissae and ordinates of a function.
  26.  *              (A minimum of four pairs is required.) It
  27.  *              produces a similar set, which is approximately
  28.  *              equally spaced and includes the input set, on the
  29.  *              standard output. The cubic spline output (R.W.
  30.  *              Hamming, "Numerical Methods for Scientists and
  31.  *              Engineers", 2nd ed. 349ff) has two continuous
  32.  *              derivatives and sufficiently many points to look
  33.  *              smooth when plotted.
  34.  *
  35.  *              The following options are recognized, each as a
  36.  *              separate argument:
  37.  *
  38.  *              -a  Supply abscissae automatically (they are
  39.  *                  missing from the input); spacing is given by
  40.  *                  the next argument or is assumed to be 1 if
  41.  *                  next argument is not a number.
  42.  *
  43.  *              -k  The constant "k" is used in the boundary
  44.  *                  value computation
  45.  *
  46.  *                      y " = ky " , y " = ky "
  47.  *                       0      1     n      n-1
  48.  *
  49.  *                  is set by the next argument. By default,
  50.  *                  k = 0. A value of k = 0.5 often results in a
  51.  *                  smoother curve at the endpoints than the
  52.  *                  default value. Negative values for k are not
  53.  *                  allowed. Cannot be used with -p option.
  54.  *
  55.  *              -n  Next argument (which must be an integer)
  56.  *                  specifies the number of intervals that are to
  57.  *                  occur between the lower and upper "x" limits.
  58.  *                  If -n option is not given, default spacing is
  59.  *                  100 intervals.
  60.  *
  61.  *              -p  Make output periodic, i.e. match derivatives
  62.  *                  at ends. First and last input values must
  63.  *                  agree. Cannot be used with -k option.
  64.  *
  65.  *              -x  Next 1 (or 2) arguments are lower (and upper)
  66.  *                  "x" limits. Normally these limits are
  67.  *                  calculated from the data. Automatic abscissae
  68.  *                  start at lower limit (default 0). If either
  69.  *                  argument is outside of the range of
  70.  *                  abscissae, it is ignored.
  71.  *
  72.  * Diagnostics: When data is not strictly monotone in "x", SPLINE
  73.  *              reproduces the input without interpolating extra
  74.  *              points.
  75.  *
  76.  * Bugs:        A limit of 1000 input points is silently
  77.  *              enforced.
  78.  *
  79.  *              The -n option has not been implemented in
  80.  *              accordance with the "UNIX Programmer's Manual"
  81.  *              specification. This was done to avoid ambiguities
  82.  *              when the -n option follows the -x option with one
  83.  *              argument.
  84.  *
  85.  *              At certain negative values for the -k option (for
  86.  *              example, k equals -4.0), the curve becomes
  87.  *              discontinuous. The -k option value has thus been
  88.  *              arbitrarily constrained to be greater than or
  89.  *              equal to zero.
  90.  *
  91.  * Credits:     The above description is a reworded and expanded
  92.  *              version of that appearing in the "UNIX Programmer's
  93.  *              Manual", copyright 1979, 1983 Bell Laboratories.
  94.  */
  95.  
  96. /*** Definitions ***/
  97.  
  98. #define FALSE        0
  99. #define TRUE         1
  100. #define MAX_SIZE  1000  /* Input point array limit */
  101.  
  102. #define ILL_ARG      0  /* Error codes */
  103. #define ILL_CMB      1
  104. #define ILL_KVL      2
  105. #define ILL_NVL      3
  106. #define ILL_OPT      4
  107. #define ILL_XVL      5
  108. #define INS_INP      6
  109. #define MIS_KVL      7
  110. #define MIS_NVL      8
  111. #define MIS_XVL      9
  112. #define MIS_YVL     10
  113. #define NMT_ORD     11
  114.  
  115. #define SQUARE(a) a*a
  116. #define CUBE(a) a*a*a
  117.  
  118. /*** Typedefs ***/
  119.  
  120. typedef int BOOL;       /* Boolean flag */
  121.  
  122. /*** Include Files ***/
  123.  
  124. #include <stdio.h>
  125. #include <ctype.h>
  126. #include <math.h>
  127.  
  128. /*** Main Body of Program ***/
  129.  
  130. int main(argc,argv)
  131. int argc;
  132. char **argv;
  133. {
  134.   int n = 0,
  135.       i,
  136.       j,
  137.       n_val = 0,
  138.       atoi();
  139.   float x[MAX_SIZE],
  140.         y[MAX_SIZE];
  141.   double a_val = 1.0,
  142.          k_val = 0.0,
  143.          x1_val = 0.0,
  144.          x2_val = 0.0,
  145.          x_intvl,
  146.          ix,
  147.          iy,
  148.          d2y[MAX_SIZE],
  149.          h,
  150.          atof(),
  151.          fabs(),
  152.          spl_int();
  153.   char buffer[257],
  154.        *temp,
  155.        *gets();
  156.   BOOL aflag = FALSE,   /* Command-line option flags */
  157.        kflag = FALSE,
  158.        pflag = FALSE,
  159.        x1flag = FALSE,
  160.        x2flag = FALSE,
  161.        is_float();
  162.   void spl_coeff(),
  163.        pspl_coeff(),
  164.        error();
  165.  
  166.   /* Parse the command line for user-selected options */
  167.  
  168.   while(--argc)
  169.   {
  170.     temp = *++argv;
  171.     if(*temp != '-')    /* Check for legal option flag */
  172.       error(ILL_OPT,*argv);
  173.     else
  174.       switch(toupper(*++temp))
  175.       {
  176.         case 'A':       /* "-a" option */
  177.           aflag = TRUE;
  178.           if(argc > 1 && is_float(*(argv+1)))
  179.           {
  180.             argc--;
  181.             argv++;
  182.             if((a_val = atof(*argv)) <= 0.00)
  183.               error(ILL_ARG,*argv);
  184.           }
  185.           break;
  186.         case 'K':       /* "-k" option */
  187.           if(pflag == TRUE)
  188.             error(ILL_CMB,NULL);
  189.           kflag = TRUE;
  190.           if(argc > 1 && is_float(*(argv+1)))
  191.           {
  192.             argc--;
  193.             argv++;
  194.             k_val = atof(*argv);
  195.             if(k_val < 0.00)
  196.               error(ILL_KVL,*argv);
  197.             break;
  198.           }
  199.           else
  200.             error(MIS_KVL,NULL);
  201.         case 'N':       /* "-n" option */
  202.           if(argc > 1)
  203.           {
  204.             argc--;
  205.             argv++;
  206.             if((n_val = atoi(*argv)) < 1)
  207.               error(ILL_NVL,*argv);
  208.             else
  209.               break;
  210.           }
  211.           else
  212.             error(MIS_NVL,NULL);
  213.         case 'P':       /* "-p" option */
  214.           if(kflag == TRUE)
  215.             error(ILL_CMB,NULL);
  216.           pflag = TRUE;
  217.           break;
  218.         case 'X':       /* "-x" option */
  219.           x1flag = TRUE;
  220.           if(argc > 1 && is_float(*(argv+1)))
  221.           {
  222.             argc--;
  223.             argv++;
  224.             x1_val = atof(*argv);
  225.           }
  226.           else
  227.             error(MIS_XVL,NULL);
  228.           if(argc > 1 && is_float(*(argv+1)))
  229.           {
  230.             x2flag = TRUE;
  231.             argc--;
  232.             argv++;
  233.             x2_val = atof(*argv);
  234.             if(x2_val <= x1_val)
  235.               error(ILL_XVL,x2_val);
  236.           }
  237.           break;
  238.         default:        /* "-n" option */
  239.           error(ILL_OPT,*argv);
  240.       }
  241.   }
  242.   if(n_val == 0)        /* Set "n_val" if not given */
  243.     n_val = 100;
  244.  
  245.   /* Get the input data */
  246.  
  247.   while(1)              /* ... while there is more input data */
  248.   {
  249.     if(aflag == TRUE)   /* Automatic abscissae were called for */
  250.     {
  251.       if(n == 0)
  252.         x[0] = x1_val;
  253.       else
  254.         x[n] = x[n-1] + a_val;
  255.     }
  256.     else                /* Abscissae supplied with input data */
  257.     {
  258.       if(gets(buffer))
  259.         x[n] = atof(buffer);
  260.       else
  261.         break;
  262.     }
  263.     if(gets(buffer))    /* Read in the corresponding ordinate */
  264.       y[n] = atof(buffer);
  265.     else
  266.     {
  267.       if(aflag == TRUE)
  268.         break;
  269.       else
  270.         error(MIS_YVL,NULL);
  271.     }
  272.     if(++n == MAX_SIZE) /* Maximum amount of input data? */
  273.       break;
  274.   }
  275.   if(n < 4)             /* Check for insufficient input data */
  276.     error(INS_INP,NULL);
  277.  
  278.   /* Check for non-monotonic abscissae. Output input data set
  279.    * without interpolation if true.
  280.    */
  281.  
  282.   h = x[1] - x[0];      
  283.   for(i = 1; i < n-1; i++)
  284.     if(fabs(x[i+1] - x[i] - h) > 0.0)
  285.     {
  286.       for(i = 0; i < n; i++)
  287.         printf("%g\n%g\n",x[i],y[i]);
  288.       exit();
  289.     }
  290.  
  291.   /* Calculate abscissa interval. Use "-x" option values if
  292.    * they were given unless they fall outside the range of
  293.    * given (or calculated) abscissae.
  294.    */ 
  295.  
  296.   if(x1flag == FALSE || x1_val < x[0])
  297.     x1_val = x[0];
  298.   if(x2flag == FALSE || x2_val > x[n-1])
  299.     x2_val = x[n-1];
  300.   x_intvl = (x2_val - x1_val)/n_val;
  301.  
  302.   /* Find the coefficients */
  303.  
  304.   if(pflag == FALSE)
  305.     spl_coeff(y,d2y,h,n,k_val);
  306.   else
  307.     pspl_coeff(y,d2y,h,n);
  308.  
  309.   /* Interpolate and output results */
  310.  
  311.   ix = x1_val;
  312.   i = 1;
  313.   for(j = 0; j <= n_val; j++)
  314.   {
  315.     while(ix >= x[i] && i < n - 1)
  316.       i++;
  317.     iy = spl_int(ix,x,y,d2y,h,i);
  318.     printf("%g\n%g\n",ix,iy);
  319.     ix += x_intvl;
  320.   }
  321. }
  322.  
  323. /*** Functions ***/
  324.  
  325. /* SPL_COEFF() - Calculate spline coefficients and return in
  326.  *               vector "d2y[]". Matrix to be solved has the
  327.  *               typical form:
  328.  *
  329.  *      +-                -+ +-   -+       +-              -+
  330.  *      |  4+k 1   0   0   | | y1" | = m * | y2 - 2*y1 + y0 |
  331.  *      |  1   4   1   0   | | y2" |       | y3 - 2*y2 + y1 |
  332.  *      |  0   1   4   0   | | y3" |       | y4 - 2*y3 + y2 |
  333.  *      |  0   0   1   4+k | | y4" |       | y5 - 2*y4 + y3 |
  334.  *      +-                -+ +-   -+       +-              -+
  335.  *
  336.  *               where k = k_val, m = 6.0/(h*h) and yn" is the
  337.  *               second derivative of the interpolated function
  338.  *               at the "nth" abscissa, y0" = k * y1" and
  339.  *               y5" = k * y4".
  340.  */
  341.  
  342. void spl_coeff(y,d2y,h,n,k_val)
  343. float y[];
  344. double d2y[],
  345.        h,
  346.        k_val;
  347. int n;
  348. {
  349.   double m,
  350.          a[MAX_SIZE-1];
  351.   int i;
  352.  
  353.   /* Set up the (symmetric tridiagonal) matrix, where the only
  354.    * elements of interest are those on the diagonal. These are
  355.    * stored in array "a[]". Array "d2y[]" initially holds the
  356.    * constants vector, then is overlaid with the calculated
  357.    * variables vector.
  358.    */
  359.  
  360.   m = 6.0/(h*h);
  361.   for(i = 1; i < n-1; i++)
  362.     d2y[i] = m * (y[i+1] - 2.0 * y[i] + y[i-1]);
  363.   a[1] = 4.0 + k_val;
  364.  
  365.   /* Reduce the matrix to upper triangular form */ 
  366.  
  367.   for(i = 2; i < n-2; i++)
  368.   {
  369.     a[i] = 4.0 - 1.0/a[i-1];
  370.     d2y[i] -= d2y[i-1]/a[i-1];
  371.   }
  372.   a[n-2] = 4.0 + k_val - 1.0/a[n-3];
  373.   d2y[n-2] -= d2y[n-3]/a[n-3];
  374.   d2y[n-2] /= a[n-2];
  375.  
  376.   /* Solve using back substitution */
  377.  
  378.   for(i = n-3; i > 0; i--)
  379.     d2y[i] = (d2y[i] - d2y[i+1])/a[i];
  380.  
  381.   /* Solve for end conditions */
  382.  
  383.   d2y[n-1] = d2y[n-2] * k_val;
  384.   d2y[0] = d2y[1] * k_val;
  385. }
  386.  
  387. /* PSPL_COEFF() - Calculate periodic spline coefficients and
  388.  *                return in vector "d2y[]". Matrix to be solved
  389.  *                has the typical form:
  390.  *
  391.  *      +-             -+ +-   -+       +-                 -+
  392.  *      | 4  1  0  0  1 | | y1" | = m * | y2 - 2*y1 + y0    |
  393.  *      | 1  4  1  0  0 | | y2" |       | y3 - 2*y2 + y1    |
  394.  *      | 0  1  4  1  0 | | y3" |       | y4 - 2*y3 + y2    |
  395.  *      | 0  0  1  4  1 | | y4" |       | y5 - 2*y4 + y3    |
  396.  *      | 1  0  0  1  4 | | y5" |       | y1 - y0 - y5 + y4 |
  397.  *      +-             -+ +-   -+       +-                 -+
  398.  *
  399.  *               where m = 6.0/(h*h) and yn" is the second
  400.  *               derivative of the interpolated function at the
  401.  *               "nth" abscissa and y0" = y5".
  402.  */
  403.  
  404. void pspl_coeff(y,d2y,h,n)
  405. float y[];
  406. double d2y[],
  407.        h;
  408. int n;
  409. {
  410.   double c,
  411.      m,
  412.      fabs();
  413.   static double a[MAX_SIZE-1],
  414.                 b[MAX_SIZE];
  415.   int i;
  416.  
  417.   /* Check for matching end ordinates */
  418.  
  419.   if(fabs(y[n-1] - y[0]) > 0.0)
  420.     error(NMT_ORD,NULL);
  421.  
  422.   /* Set up the matrix, where array "a[]" holds the diagonal
  423.    * elements, "b[]" holds those elements of column "n-1", and
  424.    * "c" holds the current element of interest of row "n-1".
  425.    * Array "d2y[]" initially holds the constants vector, then is
  426.    * overlaid with the calculated variables vector.
  427.    */
  428.  
  429.   m = 6.0/(h*h);
  430.   for(i = 1; i < n-1; i++)
  431.     d2y[i] = m * (y[i+1] - 2.0 * y[i] + y[i-1]);
  432.   d2y[n-1] = m * (y[1] - y[0] - y[n-1] + y[n-2]);
  433.   a[1] = 4.0;
  434.   b[1] = 1.0;
  435.   b[n-2] = 1.0;
  436.   b[n-1] = 4.0;
  437.   c = 1.0;
  438.  
  439.   /* Reduce the matrix to upper triangular form */ 
  440.  
  441.   for(i = 2; i < n-1; i++)
  442.   {
  443.     m = 1/a[i-1];
  444.     a[i] = 4.0 - m;
  445.     b[i] -= b[i-1] * m;
  446.     d2y[i] -= d2y[i-1] * m;
  447.     b[n-1] -= c * m * b[i-1];
  448.     d2y[n-1] -= c * m * d2y[i-1];
  449.     c = -c * m;
  450.   }
  451.   c += 1.0;
  452.   b[n-1] -= c * b[n-2]/a[n-2];
  453.   d2y[n-1] -= c * d2y[n-2]/a[n-2];
  454.  
  455.   /* Solve using back substitution */
  456.  
  457.   d2y[0] = d2y[n-1] /= b[n-1];
  458.   d2y[n-2] = (d2y[n-2] - b[n-2] * d2y[n-1])/a[n-2];
  459.   for(i = n-3; i > 0; i--)
  460.     d2y[i] = (d2y[i] - b[i] * d2y[n-1] - d2y[i+1])/a[i];
  461. }
  462.  
  463. /* SPL_INT - Interpolate points using spline function */
  464.  
  465. double spl_int(ix,x,y,d2y,h,i)
  466. float x[],
  467.       y[];
  468. double ix,
  469.        d2y[],
  470.        h;
  471. int i;
  472. {
  473.   double iy,
  474.          t1,
  475.          t2;
  476.  
  477.   t1 = (ix - x[i-1])/h;
  478.   t2 = (x[i] - ix)/h;
  479.   iy = y[i-1] * t2 + y[i] * t1 - SQUARE(h) * (d2y[i-1] * (t2 -
  480.       CUBE(t2)) + d2y[i] * (t1 - CUBE(t1)))/6.0;
  481.   return iy;
  482. }
  483.  
  484. /* IS_FLOAT() - Check that character string is in correct floating
  485.  *              point format. Return TRUE if correct, FALSE
  486.  *              otherwise. The algorithm used is a deterministic
  487.  *              finite state machine. Using the regular
  488.  *              expression terminology of Unix's "lex", the
  489.  *              character string must be of the form:
  490.  *
  491.  *                      -?d*.?d*(e|E(\+|-)?d+)?
  492.  *
  493.  *              where d = 0|1|2|3|4|5|6|7|8|9
  494.  */
  495.  
  496. BOOL is_float(str)
  497. char *str;
  498. {
  499.   int c,                /* Next FSM input character */
  500.       state = 0;        /* Current FSM state */
  501.  
  502.   while(c = *str++)
  503.   {
  504.     switch(state)
  505.     {
  506.       case 0:           /* FSM State 0 */
  507.         switch(c)
  508.         {
  509.           case '-':     
  510.             state = 1;
  511.             break;
  512.           case '.':
  513.             state = 3;
  514.             break;
  515.           default:
  516.             if(isdigit(c))
  517.               state = 2;
  518.             else
  519.               return FALSE;
  520.             break;
  521.         }
  522.         break;
  523.       case 1:           /* FSM State 1 */
  524.         switch(c)
  525.         {
  526.           case '.':
  527.             state = 2;
  528.             break;
  529.           default:
  530.             if(isdigit(c))
  531.               state = 2;
  532.             else
  533.               return FALSE;
  534.             break;
  535.         }
  536.         break;
  537.       case 2:           /* FSM State 2 */
  538.         switch(c)
  539.         {
  540.           case '.':
  541.             state = 4;
  542.             break;
  543.           case 'e':
  544.           case 'E':
  545.             state = 5;
  546.             break;
  547.           default:
  548.             if(isdigit(c))
  549.               state = 2;
  550.             else
  551.               return FALSE;
  552.             break;
  553.         }
  554.         break;
  555.       case 3:           /* FSM State 3 */
  556.         if(isdigit(c))
  557.           state = 4;
  558.         else
  559.           return FALSE;
  560.         break;
  561.       case 4:           /* FSM State 4 */
  562.         switch(c)
  563.         {
  564.           case 'e':
  565.           case 'E':
  566.             state = 5;
  567.             break;
  568.           default:
  569.             if(isdigit(c))
  570.               state = 4;
  571.             else
  572.               return FALSE;
  573.             break;
  574.         }
  575.         break;
  576.       case 5:           /* FSM State 5 */
  577.         switch(c)
  578.         {
  579.           case '+':
  580.           case '-':
  581.             state = 6;
  582.             break;
  583.           default:
  584.             if(isdigit(c))
  585.               state = 7;
  586.             else
  587.               return FALSE;
  588.             break;
  589.         }
  590.         break;
  591.       case 6:           /* FSM State 6 */
  592.         if(isdigit(c))
  593.           state = 7;
  594.         else
  595.           return FALSE;
  596.         break;
  597.       case 7:           /* FSM State 7 */
  598.         if(isdigit(c))
  599.           state = 7;
  600.         else
  601.           return FALSE;
  602.         break;
  603.     }
  604.   }
  605.   return TRUE;
  606. }
  607.  
  608. /* ERROR() - Error reporting. Returns an exit status of 2 to the
  609.  *           parent process.
  610.  */
  611.  
  612. void error(n,str)
  613. int n;
  614. char *str;
  615. {
  616.   fprintf(stderr,"\007\n*** ERROR - ");
  617.   switch(n)
  618.   {
  619.     case ILL_ARG:
  620.       fprintf(stderr,"Argument must be greater than zero: %s",
  621.           str);
  622.       break;
  623.     case ILL_CMB:
  624.       fprintf(stderr,"Cannot use -k option with -p option");
  625.       break;
  626.     case ILL_KVL:
  627.       fprintf(stderr,"Illegal value for -k option: %s",str);
  628.       break;
  629.     case ILL_NVL:
  630.       fprintf(stderr,"Illegal value for -n option: %s",str);
  631.       break;
  632.     case ILL_OPT:
  633.       fprintf(stderr,"Illegal command line option: %s",str);
  634.       break; 
  635.     case ILL_XVL:
  636.       fprintf(stderr,"Illegal value for -x option: %s",str);
  637.       break;
  638.     case INS_INP:
  639.       fprintf(stderr,"Insufficient input data");
  640.       break; 
  641.     case MIS_KVL:
  642.       fprintf(stderr,"Missing value for -k option");
  643.       break;
  644.     case MIS_NVL:
  645.       fprintf(stderr,"Missing value for -n option");
  646.       break;
  647.     case MIS_XVL:
  648.       fprintf(stderr,"Missing value for -x option");
  649.       break;
  650.     case MIS_YVL:
  651.       fprintf(stderr,"Missing ordinate value");
  652.       break;
  653.     case NMT_ORD:
  654.       fprintf(stderr,"End ordinates do not match");
  655.       break;
  656.     default:
  657.       break;
  658.   }
  659.   fprintf(stderr," ***\n\nUsage: spline [-aknpx]\n");
  660.   exit(2);
  661. }
  662.  
  663. /*** End of SPLINE.C ***/
  664.